% 根据模型长度和密集度，生成质量阵、刚度阵
% [nBall] ~o~~o~~o~~o~
% [DOF]  .~o~~o~~o~~o~

function [M,K,C,m,k,nBall,DOF]=func_model(L,ctn,cm,ck)

fprintf('[ODC] func_model: 正在生成 K, M, C ...\n')

nBall=L*ctn; 
DOF=nBall+1; % 显式计算，边界条件施加方式。需要多出一个点"."，作为激励作用点。

m=1/ctn;
k=1*ctn;
l0=L/nBall; %小球间隔，实际计算中用不到

% 质量阵（稀疏存储）
% 注意此处，激励作用点也赋予了质量m，这在能量计算中会多出一部分，记得去掉。
vector_m=m*ones(DOF,1);
M=spdiags(vector_m,0,DOF,DOF);

% 刚度阵（稀疏存储）
vector_k=k*ones(DOF,1);
K=spdiags([-vector_k,2*vector_k,-vector_k],[-1,0,1],DOF,DOF);
K(1,1)=2*k;
K(1,2)=-2*k;
K(2,1)=-2*k;
K(2,2)=3*k;
K(DOF,DOF)=k;

% 生成阻尼阵. 此处存疑：或是? C=cm*ctn*M+ck/ctn*K;
C=cm*M+ck*K;

end